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We numerically generate, and then study the basic properties of dark soliton-like excitations in 
a dipolar gas confined in a quasi one dimensional trap. These excitations, although very similar 
to dark solitons in a gas with contact interaction, interact with each other and can form bound 
states. During collisions these dipolar solitons emit phonons, loosing energy, but accelerating. Even 
after thousands of subsequent collisions they survive as gray solitons and finally reach dynamical 
equilibrium with background quasiparticles. Finally, in the frame of classical field approximation, 
we verified, that these solitons appear spontaneously in thermal samples, analogously to the type II 
excitations in a gas of atoms with contact interaction. 


PACS numbers: 03.75.Hh , 03.75.Lm , 05.45.Yv 

I. INTRODUCTION 

The main motivation to our work comes from re¬ 
cent progress concerning the role of dark solitons in the 
physics of the gas with contact interaction only. First, 
let us remind, that solitons are localized waves propa¬ 
gating through non-linear media without any spreading 
and also moving through each other without changing 
the shape. Among equations of broad interest for physi¬ 
cists, there are a few which exhibit solitons. One of them 
is the non-linear one dimensional Schrodinger equation: 

ihdt ^{z) = (1) 

in which one can find bright (the case with attracting 
forces, ^ < 0) and dark (the repulsive case ^ > 0) solitons. 

The Eq. ([1]) applies, among others, to the dynamics of 
bosons interacting via a (i-like potential, i.e. the interac¬ 
tion potential between two bosons V{z — z') is equal to 
g6{z — z'). In this context, the Eq. ^ is only an approx¬ 
imation valid in the variously defined weak interaction 
limit [ij. Such a gas of bosons, can be described by more 
fundamental many-body model, the Lieb-Liniger model, 
which has been used to find rigorously the ground state 
(H as well as the elementary excitations It turns out, 
that there are two branches of elementary excitations, 
staying on the same footing. The first branch has been 
identified as the Bogoliubov excitations. The second one, 
consisting of the so called ’’type II” excitations, within 
research performed by the next generations of physicists 
turned out to be solitonic branch M- It is now clear 
that the dark solitons, as the elementary excitations, ap- 
pear spontaneously in a gas at finite temperature (1Q| |. 

The equation dm, is not capturing many important 
and common situations, even in the weakly interacting 
limit. For example, it is neglecting possible external 
traps in which the gas is confined, and it is assuming 
contact interaction only. In particular, it misses the gas 
with a dipolar interaction, which is recently studied in 
more and more laboratories [nHi3. Among many moti¬ 
vations to study the dipolar gas is definitely its unusual 


excitation spectrum, which in certain situations, close to 
instabilities, can have local minimum for non-zero quasi 
momenta, so called roton spectrum (a. It seems natural 
to ask if here is a place for the second branch and from 
what excitations it can be build. Up to our knowledge, 
there is no study of the elementary excitation spectrum 
beyond the Bogoliubov approximation, neither by solv¬ 
ing the underlying many body problem, nor on the mean 
field level. 


Moreover, although the literature about the bright soli¬ 
tons is very reach (see parts of the review 0), relatively 
little is known about dark solitons. For the time being 
the shapes of dipolar solitons in a gas with both con¬ 
tact and dipolar interaction has been presented [il0. 
Surprisingly, the life time of the dark solitons can be in¬ 
creased with some tricks even in the higher dimensional 
space [1^. The dark solitons in non-local nonlinear media 
were the subject of research in optics. It has been found 
that two optical dark solitons interact with each other 
via non-trivial potential, they can even form a bound 
state (lOl. In this context, although the effective forces 
between photons can be non-local, it is always assumed 
that they are not long-range. 


Within this work, limited to the mean field descrip¬ 
tion, we look closer at the basic properties of the dark 
solitons in the dipolar media. The aim of this paper is 
to study dark dipolar solitons and contrast them with 
the true solitons existing in the gas with contact inter¬ 
action only. The analysis is performed within the mean 
field description, but we include also the finite tempera¬ 
ture case. We start with the detailed description of our 
model in Sec. nn After presentation of the basic features 
of single soliton (shapes, dynamics) in Sec IIIIl we discuss 
a motion of two dipolar solitons in Sec. IIVI Sec. [Vl is 
devoted to the analysis of the thermal samples. 
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Figure 1: (Color online) Configuration: Trapping potential is 
given by a steep harmonic trap in the X and the Y directions 
(called transverse directions), and a long box with periodic 
boundary conditions in the Z direction. The dipole moments 
of the constituent atoms are polarized along X axis, such that 
they mostly repel each other. 


II. THE MODEL 

We assume that the dipolar gas is confined in the po¬ 
tential 


U{x,y,z) = ]^rnw\{x^+ y^), (2) 

with a tight harmonic trap in the X and the Y directions. 
The space in the Z direction is assumed to be finite, 
with the length L and with periodic boundary conditions. 
These conditions impose a quantization of the momenta 
kz in this direction, kz is a multiple of 

Throughout the paper we use the box units, where L, 
ml? k? jml? and 1?/ml?k^ are the units of length, 
time, energy and temperature, respectively. We assume 
that the dipole moments of the constituent atoms are 
polarized along the X axis , namely when the dipolar 
potential in the momentum space has the form jsH 


Kiip = 


^dd 


(2^) 


3/2 



(3) 


where is a dipolar coupling strength. Thus, the 
atoms move on the circumference of a circle, having the 
dipole moments perpendicular to the circle-plane. In this 
configuration dipoles mostly repel each other, what fa¬ 
vors the dark solitons. The geometry of our problem is 
sketched in Fig. [H 

We will investigate the system within the mean field 
model. The time-dependent mean field wave-function 
obeys the three dimensional Gross-Pitaevskii equation 
(GPE) 


ihdt'ipw = {-- - VU{x,y,z)+gN\'ip:iu?‘ 

2m 

+ J dr'\^3B{r')\‘^Vdip{r - r'))'ip3B{r). 

We assume that the transverse confinement is so tight, 
that the wave-function stays in the lowest energy level 


in the X and Y directions, namely when the following 
single-mode approximation holds: 

'>p 3 D{x,y,z;t) = e~''‘^^^(f>o{x)(j)o{y)ip{z,t), (4) 

where is the Gaussian wave function: 



where l± = y/h/{muj±) is the transverse oscillator width. 
The necessary condition underlying this approximation is 
that both, the thermal energy /cpT and the interaction 
energy, are much below the energy of the first excited 
state in the transverse direction. 


ksT^gn huj±. ( 6 ) 

Assuming such a situation we multiply the three dimen¬ 
sional GPE by (j)o{x)(f)o{y), integrate the result with re¬ 
spect to X and y and finally we reach the one dimensional 
effective GPE 

indti! = ^ j dz'Vee{z') ipiz) (7) 

The Eourier transform of the effective one dimensional 
dipolar potential Ve^ik) reads [illlil. 

VMk) = g,on + 3ffdd (1 + / {{XkLf/2)), (8) 

where the constant term with ^con = ^ mimics a 

contact interaction, we introduced one dimensional dipo- 

3D 

lar coupling strength and X = l±/L is equal 

to the aspect ratio of the trap. The function / which 
appears in Eq. (|8|) is equal to /(a) = ae"Ei(a), where 
Ei is the exponential integral. In this paper we aim at 
studying the role of pure dipolar interactions only, hence 
from now on, we will assume ^con = 0 (d^ . 

The spatial representation of this potential is visual¬ 
ized in Eig. [21 It is worth stressing that Eq. © describes 
not a movement of single dipoles, but rather the move¬ 
ment of slices of the dipolar gas as marked in Eigs. (Hand 
m If such slices are far from each other, then the details 
of the distribution of dipoles within such slices plays no 
role, and we can treat them as single large dipoles. In 
such situations the interaction between them scales as 
1/z^. The situation changes, when the distance between 
them is of the order of the transverse width lj_. This gives 
the characteristic range of the effective ID potential ([8]). 
Einally in the limit of slices going to the same place, the 
effective dipolar potential converges to a constant. 


III. SINGLE SOLITONS 

In the gas with contact interaction only, the phase of 
a solitonic wave function is changing by tt around the 
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Figure 2: (Color online) Effective potential Ves defined in (|8|) 
in position representation. For large distances z it behaves 
like 1/z^, its characteristic range (width at half maximum) is 
of the order of l±. Parameters: ^con = 0 and gad = 1- 



Figure 3: (Color online) Shape of the solitonic-like excitations 
in the dipolar gas. Results of the imaginary time evolution 
with the phase constraint given in Eq. are marked with: 
red solid line for A = 10“^, dashed blue for A = 10“^ and dot- 
dashed green for A = 10“^. fn all figures: Ngaa = 4000 and 
gcon = 0. Gray thick line indicates the Shabat-Zakharov dark 
soliton for the gas with contact interaction only, Ngcon = 
12000. The blue dashed line corresponds to experimental-like 
parameters: the mass and the dipole moment of Dysprosium, 
assuming N = 600 atoms and the transverse trap equal to 
uj±_ = 27r X 5kHz. 


region with minimal density. This motivates us to look 
for the solitons in the dipolar gas within the wave func¬ 
tion which minimizes the energy but under constraint of 
a jump of the phase equal to tt. If we only force the jump 
of the phase equal to tt, then the periodic boundary con¬ 
dition would not be satisfied. Hence, except the tt jump 
required by the soliton, we add to the constraint on the 


phase a linearly changing function of the position: 

= \{^- sign(2:)^ , (9) 

such that the phases at the edges of the box are equal. 
Due to this additional term the gas is moving with the 
speed tt/L. When showing the result for the dynamics, 
we would unrotate this motion. 

To find the state with the lowest energy but with the 
phase equal to the constraint (|9]) we use the method of the 
imaginary time evolution, forcing at each step the phase 
®. Densities computed with this method for three dif¬ 
ferent aspect ratios A = 10“^, 10“^ and 10“^ are shown 
in Fig. [3l We compare these profiles with the density of 
a dark soliton from purely contact interacting case (22| 
marked in Fig[3]with the thick gray line. As in both cases, 
a contact interacting gas and the considered here quasi ID 
dipolar gas, the lowest lying excitations are phonons, 
then comparing their speed of sound we translate the 
dipolar coupling coefficient ^fdd into an effective coupling 
coefficient ^eff- We find ^eff = ^9dd- The speed of sound 
in the dipolar media is equal to \^3rigddJ^- The width of 
the Shabat-Zakharov soliton is proportional to the heal- 
ing length, f 

where Ueff is the effective scattering length defined by 
the relation ^eff = The Gross-Pitaevskii equa¬ 

tion is supposed to be valid in the weakly interacting 
limit, namely when nueff ^ 1, what implies ^ I±. 

One can see that the more elongated is the trap, the 
closer is the result to the soliton in contact interacting 
gas, marked in Fig. [3] with the thick gray line. To un¬ 
derstand this observation mathematically we look at the 
system in the momentum space. The range of the effec¬ 
tive dipolar potential in the momentum space is equal 
to 1//_L, equal in the box units to 1/A. Hence, in the 
limit of a very elongated trap, A 0 the range of the 
effective potential is going to oo, and the evolution is 
dominated by its values relatively close to k = 0. The 
derivative of the effective potential Eq. m ai k = 0 van¬ 
ishes, what means that the potential is flat at the origin. 
Thus locally, where this ’local’ region spans even to high 
momenta, the interaction potential looks like a constant 
function, as it would be in the contact interacting case. 
The ratio between the typical width of the soliton and the 
range of the potential , written in the momentum space 
and in the box units is equal to ^ = ^/ng^X. Thus if the 
trap is very elongated and the coupling strength is small, 
we expect that the dipolar soliton will be close to the 
Shabat-Zakharov solution and moreover that our system 
will evolve similarly to the gas with contact interactions 
only. The presented solutions are only candidates for 
solitons, having similar jump of the phase as the Shabat- 
Zakharov solitons. The crucial feature of solitons is their 
ability to propagate without changing their shape. To 
verify if this is the case also here, we were setting as 
initial states the results of the imaginary time evolution, 
and trace their dynamics obtained via numerical solution 
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Figure 4: Density as a function of time presented in the frame 
rotating with the speed tt/L. In the laboratory frame the 
solitons travel around the circle 17 times. Parameters: A = 
0.01, ngdd = 4000 and n^con = 0, as in the Dysprosium like 
case commented in caption of Fig. O The horizontal red line 
is at time tphon = 0.16 : a time at which a single phonon 
would pass the whole box. 


of GPE. An example is given in Fig. IH where we unro¬ 
tated the constant speed due to the linear dependance of 
the phase. The evolution time in this figure is about 30 
times longer than a time tphon at which a phonon would 
travel across the whole box. When testing the system in 
a wide parameter range and for a longer evolution time 
we always observe the dips which propagate without not- 
icable distortion. 


IV. TWO SOLITONS 



POSITION (box units) 

Figure 5: Density as a function of time. The initial distance 
between solitons is equal to 0.4. Parameters as in FigO 

The real solitons obey an equivalent of the superposi¬ 
tion principle: although the system is non-linear one can 


construct many solitonic solution based on a single one. 
As shown in the paper of Shabat and Zakharov (22| in 
the contact interacting gas two solitons, when colliding, 
do not change the shape. To test if the dipolar solitons 
can behave similarly, we generate and propagate pairs of 
them. We follow the method described for a single soli- 
ton, but this time we enforce two jumps of the phase, at 
positions 2:1 and 2 : 2 , and without linear background 

TT 

^ 2 {z) = - (sign( 2 : - zi) - sign( 2 ; - Z 2 )). (10) 

Having such constraint on the phase we reach with the 
help of the imaginary time evolution two solitons and 
then we propagate them using GPE. An example of such 
dynamics is shown in EigO The immediate conclusion is 
that the dipolar solitons interact , in the example shown 
they attract each other. Hence, the excitations we study 
are not solitons in the strict mathematical sense. Even 
though, these excitations resemble the shape for a long 
evolution times, thus we will use the name soliton-like or, 
for brevity, even the name soliton. 



z/L 


Figure 6: Acceleration versus distance between solitons 
asoii^Zsoi) for different initial states. Parameters as in Fig 

El 


To understand the motion of such a pair we monitor 
the position of one of the solitons versus time: Zso\{t). 
Then we numerically compute both the velocity and the 

dz^ 

acceleration asoi(t) = Combining these quanti¬ 

ties we obtain acceleration versus inter-solitonic distance 
aso\{^Zsoi), where Azsoi = 2|2:soi|- We repeat this proce¬ 
dure for many initial states. One can see that results of 
many simulations are lying on a common curve, see Fig. 
[6l It suggests us to use the following classical mechan¬ 
ics picture, in analogy to many cases for dark solitons in 
the contact interacting gas (^, treat each soliton 

as a point-like (negative) mass obeying the Newton law. 
Assuming a^oi = Egoi/m = — we can numerically 

integrate Usoi to extract inter-solitonic potential divided 
by a mass Vso\/m. We show the results for two cases: 
in Fig. [3 for very elongated trap with A = 10“^ and 
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Figure 7: Effective interaction potential between solitons Vsoi 
divided by mass in the classical picture. The figures in the 
bottom panel show density versus time for two initial con¬ 
ditions, for the minimum of the potential Azsoi{t = 0) = 
0.1232L (left) and for Azso\{t = 0) = 0.35L (right). Parame¬ 
ters: A = 0.001, Ngdd = 1000 and gcon = 0. 
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Figure 8: Effective interaction potential between solitons Esoi 
divided by mass in the classical picture. The potential con¬ 
sist of repulsive exponentially decaying core and a tail that 
decays like 1/z^. The figures in the bottom panel show den¬ 
sity versus time for three initial conditions, at the repul¬ 
sive core Azso\{t = 0) = 0.022L (left), close to the mini¬ 
mum of the potential Azso\{t = 0) = 0.031L (center) and for 
Azso\{t = 0) = 0.3L (right). Parameters as in Fig[3l 


weak dipolar interaction A^^dd = 1000 and in Fig. [8] 
for A = 10“^ and = 4000. In both cases the inter- 
solitonic potential has a repulsive core and the tail scaling 
with the distance as 1 / 2 :^. We tested this classical picture 
starting again from different initial states and watching 
their evolution. Samples of this investigation are given at 
the bottom panels in Figs. [7|and[8l Two solitons placed 
close to each other strongly repel each other. As the 
inter-solitonic potential has a local minimum, thus there 
is a possibility for bound states as shown in the panels 
marked with the letter A. Two far away lying solitons, 
like in the panel marked with the letter B, they hardly 
see each other, but after longer time they finally start 
oscillations in the common potential. 

In Figs. [3 and [8] we compare the resulting potential 
Vso \/tu just with the rescaled density of a single soliton 
a — po), where po is the density far from the 

soliton, approximately equal to 1/L (s^l- The agreement 
between solitonic shape and inter-solitonic potential is 
surprisingly good. This indicates the origin of the inter- 
solitonic potential: each soliton is simply moving in the 
mean field potential created by its partner. As the dark 
solitons have negative masses, they are climbing towards 
the local maxima of the density. 

The attentive Reader could notice some peculiarity in 
Fig. [9] bottom right panel. There, the solitons seem to 
slightly accelerate. This counter-intuitive effect is not a 
numerical artefact. To illustrate clearly what is happen- 
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Figure 9: Density as a function of time. The red solid line 
indicates the speed of sound. Parameters: A = 0.1, ngdd = 
1000 and ng con — 0. 


ing we present the evolution for the extreme conditions, 
A = 0.1 ngdd = 1000, in Fig. [9l Then one can see 
that the collision is not elastic: colliding solitons radiate 
phonons. In consequence the energy of solitons decays, 
and the dark solitons are transformed into the gray soli¬ 
tons. As in the case with contact interacting gas, the 
gray solitons have some characteristic speed. This extra 
velocity is thus gained by the solitons due to the inelastic 
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collision. 

Hence, the dissipation during the collision leads in fact 
to acceleration of solitons, as clearly shown in Fig. [9l Ac¬ 
tually this effect, although not directly described, seems 
to exist for the optical solitons with nonlocal nonlinearity 
studied in (23| (Fig. 3 therein). Until collision the move- 



TIME (box units) 


Figure 10: The grayness h defined in the text as a function 
of time. The timescale is so long, that during evolution there 
occurs hundreds of the collisions between solitons. The super¬ 
imposed red line corresponds to travel at the speed of sound. 
The plot of density vs time corresponds to the evolution at 
short time window, from t = 40 to t = 40.01. The red line is a 
fit in which we invert the energy to grayness H relation H(E) 
and assume energy E — Eq — Asinh(Ht) with fitting param¬ 
eters Fo, A and B, to model the process of the stimulated 
emission. Parameters as in Fig[3l 

ment of a pair of solitons can be captured by a classical 
model. Then, collision goes clearly beyond the classi¬ 
cal picture: it is an event when the part of the energy 
is transferred from solitons into initially unoccupied ex¬ 
citations, i.e. phonons. What is the fate of these soli¬ 
tons? To find the answer we let the system evolve for a 
long time. During the evolution we monitored minimal 
density in the system, called here a grayness H. In the 
case of a gas with contact interactions only, this param¬ 
eters can be translated into the Shabat-Zakharov soli- 
ton speed, equal to v{H) = and its energy, equal 

to E{H) = 0. The evolution is 

shown in Fig. ifol As we started with two black soli- 
tons, initially H = 0. Then, at short times,the grayness 
h increases, first quadratically and then exponentially. 
At a very long time the quantity h starts to saturate. 
We understand the dynamics in terms of interaction be¬ 
tween the subsystem (solitons) with an environment (all 
phonons and other quasiparticles), taking the idea from 
the paper (25|. Initially the environment is empty: this 
is the regime of spontaneous emission, where a coupling 
between subsystem and the unoccupied modes leads to 
a slow dissipation. Then, we enter the regime of stimu¬ 
lated emission, where the modes after the stage of spon¬ 


taneous emission are already occupied, and hence the 
rate of phonon emission rapidly grows. If the environ¬ 
ment would be infinitely large then eventually we would 
loose the energy from the subsystem completely. How¬ 
ever here, we deal with a small environment, in which the 
phonons always overlap with solitons. This allows for the 
dynamical equilibration, in which the amount of the en¬ 
ergy radiated to phonons is equal to the energy absorbed, 
as in the qubit interacting with a thermal reservoir. In¬ 
deed we see in the inset of Fig. [TOl that even after a 
long evolution time the two solitons survive. The large 
fluctuations of H are due to overlaps between phonons 
and solitons, but also due to statistical fluctuation of the 
flow of the energy between solitons and phonos. 



Figure 11: Density as a function of time at three different 
temperatures. The thick black line depicts a travel with the 
speed of sound. Parameters: A = 0.01, ng^d = 4000 and 
ngcon — 0 and N — 1000. 


V. THERMAL SOLITONS 

As mentioned in the introduction, the ultracold gas 
with contact interaction possesses not only phonons but 
also dark solitons [a. These solitons are the elementary 
excitations forming the second branch of excitations de¬ 
rived in 1963 by Elliot Lieb Q. We would like to check 
if non-Bogoliubov excitations exist also in a quasi ID 
dipolar gas. To fully verify it, one should solve a many- 
body problem, an equivalent of the Lieb-Liniger model 
for a dipolar gas. This task is definitely beyond this pa- 
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per. Instead, still relying on the mean field model, we 
can at least check if the solitons, as candidates for ’’type 
II” excitations in a dipolar gas, appear spontaneously at 
finite temperatures and if their number obeys the ther¬ 
mal statistics. At the first glance, the presence of such 
solitons seem unlikely: the system is not integrable and 
we have already shown that the dipolar soliton-like exci¬ 
tations are not robust, as they collide inelastically. On 
the other hand, we also demonstrated that two solitons 
immersed in the gas full of phonons can stay in the dy¬ 
namical equilibrium. The latter observation encourages 
us to look also at the thermal equilibrium. 

As in a, we use the classical field approximation 
(CFA), in which one replaces the quantum field opera¬ 
tor T(z) with the classical field T(z). This method is 
commonly used at zero-temperature, where ^( 2 :) is in¬ 
terpreted as the macroscopically occupied (by all parti¬ 
cles) single particle orbital, called the GPE wave func¬ 
tion. The idea of using the same replacement to describe 
a gas at non-zero temperature relies on the fact that for 
the degenerate gas there are many substantially occupied 
single-particle orbitals, which sum up to a classical field 

To compute the values of target variables averaged over 
a statistical ensemble one performs averaging over the ap¬ 
propriately chosen set of the classical fields or/and 
time-aver aging over the fields propagated with the help 
of the GPE equation or its generalization. 

There is variety of different classical field methods de¬ 
vised to describe ultracold gases (26| . Here we use a sim¬ 
ple version described in [^. Eirst, we generate many 
initial states using a method based on the Bogoliubov 
description. As our method of generating thermal sam¬ 
ples is not sufficiently good, we evolve the field with the 
help of GPE, we trace the system for sufficiently long 
time and compute time averages. 

More precisely, in each realization, the initial wave- 
function reads 

= E (11) 

^max ^max 

where the wave vector is k = 2^1 jL with I an integer. 
The sum over wave-vectors is limited by the /cmax = 
yJck^T + jjL with c = 0.29, where T is the temperature 
and /i = = 0) is the chemical potential. The 

cut-off /cmax IS choscu such that with the classical field 
method one would recover the quantum results at least 
for the ideal gas [28|, [ 2 ^ . The expansion coefficients ak 
should be drawn from the thermal classical distribution. 
Even on the classical level this drawing can be a cum¬ 
bersome task. The details of our procedure of drawing 
initial set of au are given in the Appendix. As the pro¬ 
cedure gives unsatisfactory effect (especially for higher 
temperatures), we have to propagate each guessed clas¬ 
sical field over a long time to allow for thermalization. 
Einally, after this auxiliary stage of thermalization, in 
the next stage we further propagate the state to perform 


averaging over time. The short-time windows of the evo¬ 
lution during the last stages are illustrated in Eig[TTl We 
see there three examples of the evolution of the density 
of the gas, with the fraction of the 0-momentum mode 
from 0.68 (top panel), through 0.17 (middle panel), 0.08 
(bottom panel). 

The speed of sound is depicted with the solid black 
line. In the middle and the bottom panels one can see 
dips in the density, propagating much slower than the 
speed of sound. We checked that these dips were stable 
even at long timescales. Thus we expect that they are 
dark solitons. 



< I ttg \^> 

Figure 12: Number of dark solitons versus fraction of atoms 
in the 0-momentum mode. Histogram corresponds to the 
numerically treated dipolar case with parameters A = 0.01, 
ngdd — 4000 and ngcon = 0. The green points are computed 
analogously to the dipolar case, but for ngcon — 12000 and 
gdd = 0. Finally the blue dashed line is the expected number 
of solitons assuming the equipartition of energy and assum¬ 
ing that dark solitons as quasi-particles, namely number of 
solitons is approximated with /cbT/EsoI, where T is the esti¬ 
mated temperature and Egoi is the estimated energy of the 
dark soliton. The details of the estimation are in the text. In 
each case we assumed N — 1000 atoms. The results for the 
dipolar gas and the gas with contact interaction comes from 
time averaging performed on 1000 realizations of the classical 
field, each drawn as described in the Appendix. 

To convince ourselves we prepare the following test: 
for each thermal sample we performed a long evolution, 
during which we average over time the fraction of the 
least energetic 0-momentum mode (|aoP/A^), the fluc¬ 
tuation of this fraction A (|aoP/A^) and the number of 
the deepest dips, with the density below 0.05/L. The 
last quantity, called the number of dark solitons, is pre¬ 
sented in Fig. [T 2 I We estimated the temperature of each 
sample using formulas for average number of motionless 
particles, i.e. (|(aop/A/') , for the ideal gas [sOl- As the 
crosscheck, we compute what would be fluctuation of the 
latter average in the ideal gas and compare it with the 
previously numerically found A {\ao\‘^/N). The agree¬ 
ment was within a few percent, hence we concluded that: 
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1) the interactions in the gas are still so small that the 
statistics seems close to the ideal gas case and 2) our 
estimation for temperature is reasonable within a few 
percent. We also estimated the energy of dark solitons, 
relying on the formulas for the gas with contact interac¬ 
tions only 0 as done in the previous sections. This gives 
the estimation for the energy of the dipolar dark soliton 
^soi = (ia box units). Finally, having the 

estimated energy of the dark soliton, and estimated tem¬ 
perature of each sample we computed the estimated num¬ 
ber of dark solitons /cbT/£’soi, assuming their equiparti- 
tion as it should be for quasiparticles. This expected 
number of solitons is given by the blue dashed line in 
Fig. [121 We also repeat all the steps for the gas with the 
contact interaction only. 


VI. CONCLUSIONS 

We investigated the dipolar gas confined in the po¬ 
tential having the shape of a long tube, with the dipole 
moments polarized perpendicular to the long axis. In 
this geometry we find robust localized excitations, very 
similar to the Shabat - Zakharov solitons appearing in 
the gas with contact interaction only. In contrast to the 
contact interacting case, pairs of these soliton-like dipo¬ 
lar excitations mutually interact, and are able to form 
two-soliton bound states. Moreover, the collisions be¬ 
tween such pairs are not elastic: during collision they 
emit phonons, hence they loose energy and because of 
this they accelerate. After thousands of subsequent col¬ 
lisions they can reach dynamical equilibrium with the al¬ 
ready emited quasiparticles. Finally we show that these 
excitation appear spontaneously at finite temperatures 
and their number obeys thermal statistics. The latter 
observation is closely related to the contact interacting 
case, where the appearance of the dark solitons at the 
non-zero temperatures is an evidence that they are the 
non-Bogoliubov elementary excitations, predicted by El¬ 
liot Lieb in 1963. 

There are many open questions concerning these ex¬ 
citations. One should check their stability in the real 
three dimensional calculation. It is interesting if they ex¬ 
ist also in the roton regime. There is the question of the 
effect of a harmonic trap along the long axis (here the 
box with periodic boundary was assumed). One needs 
also benchmarks from other methods used for the finite- 
temperature simulations, as they are quite developed for 
the dipolar gases (HHaii. Probably the most important 
task would be a verification if these excitations are re¬ 
ally elementary and if they are present in the many-body 
model. 
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VII. APPENDIX: GENERATING THERMAL 
STATES WITHIN CLASSICAL FIELD 
APPROXIMATION 


In the Bogoliubov method one approximates the grand 
canonical Hamiltonian with H—fiN = Ekh\hk^ where 


Ek 



^+2nyff 


is the Bogoliubov spec¬ 


trum and the operators bk and 6^. are the annihilation 
and creation bosonic ladder operators for quasiparticles, 
respectively. The operator bk, annihilating a quasiparti¬ 
cle with quasi-momenta k is related to the bosonic oper¬ 
ators d/c annihilating a particle with a given momenta k 
via relation: 


bk = dkUk - dlvk, 


( 12 ) 


where Uk = ~ 

_ 1^ classical version, we 

replace dk, d|., bk and 6^. with the complex num¬ 
bers denoted as ak, Pk ^.nd respectively. 

We introduce the classical probability distribution 

= iexp {-{T.-k^^^<k<k^.^Pk\Pk?)) 

with the partition function defined such that 
1 = j (f ... j (fak^^^P. Within this 

approximation the number of quasiparticles with mo¬ 
menta k is given by the exponential distribution with 
the parameter k^TjEk. To generate the initial state foe 
a given temperature and for a given number of atoms N 
we follow the steps: 


1) For /c 7 ^ 0 we draw the number of quasiparticles 
^/c = |/^/cP using the exponential distribution. 

2) For /c 7 ^ 0 we draw a phase (pk with the uniform 
distribution on the interval [—7r,7r). 

3) We construct the amplitude pk = y/rike^^^. 

4) We compute the amplitudes ak = Ukh E "^kP-k' 

5) We compute the amplitude of the 0-momentum 
mode, ao = 

6) We construct the classical field according to Eq. 

(na). 
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